function dmrc = dmr_CoordinateFrame(hfile, opts)
% DMR::CoordinateFrame  - generate coordinate matrices of DMR
%
% FORMAT:       dmrc = dmr.CoordinateFrame([opts]);
%
% Input fields:
%
%       opts        optional struct with settings
%        .dirs      compute directions (default true)
%        .first     compute coordinate of first voxel (default false)
%        .ndgrid    run ndgrid with DMR resolution (default false)
%        .origin    compute origin (center of STC data, default true)
%        .trans     compute transformation matrix (default true)
%        .transinv  compute inverse matrix (default false)

% Version:  v0.7b
% Build:    7083009
% Date:     Aug-30 2007, 9:32 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'dmr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'dirs') || ...
   ~islogical(opts.dirs) || ...
    isempty(opts.dirs)
    opts.dirs = true;
else
    opts.dirs = opts.dirs(1);
end
if ~isfield(opts, 'first') || ...
   ~islogical(opts.first) || ...
    isempty(opts.first)
    opts.first = false;
else
    opts.first = opts.first(1);
end
if ~isfield(opts, 'ndgrid') || ...
   ~islogical(opts.ndgrid) || ...
    isempty(opts.ndgrid)
    opts.ndgrid = false;
else
    opts.ndgrid = opts.ndgrid(1);
end
if ~isfield(opts, 'origin') || ...
   ~islogical(opts.origin) || ...
    isempty(opts.origin)
    opts.origin = true;
else
    opts.origin = opts.origin(1);
end
if ~isfield(opts, 'trans') || ...
   ~islogical(opts.trans) || ...
    isempty(opts.trans)
    opts.trans = true;
else
    opts.trans = opts.trans(1);
end
if ~isfield(opts, 'transinv') || ...
   ~islogical(opts.transinv) || ...
    isempty(opts.transinv)
    opts.transinv = false;
else
    opts.transinv = opts.transinv(1);
end

% get resolution
rsx = bc.InplaneResolutionX;
rsy = bc.InplaneResolutionY;
rsz = bc.SliceThickness + bc.GapThickness;
dmx = bc.ResolutionX;
dmy = bc.ResolutionY;
dmz = bc.NrOfSlices;

% get slice centers
cs1 = [bc.Slice1CenterX; bc.Slice1CenterY; bc.Slice1CenterZ];
csn = [bc.SliceNCenterX; bc.SliceNCenterY; bc.SliceNCenterZ];
css = mean([cs1, csn], 2);
drs = csn - cs1;
drs = drs ./ sqrt(sum(drs .* drs));

% get directions
drx = [bc.RowDirX; bc.RowDirY; bc.RowDirZ];
dry = [bc.ColDirX; bc.ColDirY; bc.ColDirZ];
drx = drx ./ sqrt(sum(drx .* drx));
dry = dry ./ sqrt(sum(dry .* dry));
drz = cross(drx, dry);
if sum(drs .* drz) < 0
    drz = -drz;
end
if abs(sum(drs .* drz) - 1) > 0.01
    error( ...
        'BVQXfile:InvalidObject', ...
        'Fields Slice*Center* and ***Dir* don''t match.' ...
    );
end
if abs(sum(drs .* ((csn - cs1) ./ (rsz * (dmz - 1)))) - 1) > 0.01
    error( ...
        'BVQXfile:InvalidObject', ...
        'Fields NrOfSlices, Slice/GapThickness, and directions mismatch.' ...
    );
end
trf = [[drx, dry, drz, css]; [0, 0, 0, 1]];

% build output
dmrc = struct;
dmrc.DimX = dmx;
dmrc.DimY = dmy;
dmrc.DimZ = dmz;

% set directions
if opts.dirs
    dmrc.DirX = drx;
    dmrc.DirY = dry;
    dmrc.DirZ = drz;
end

% compute origin
if opts.first
    fst = cs1 - ((rsx - 1) * dmx / 2) * drx + ((rsy - 1) * dmy / 2) * dry;
    dmrc.First = fst;
end

% ndgrid
if opts.ndgrid
    [dmrc.NDgrid{1:3}] = ndgrid(0:(rsx - 1), 0:(rsy - 1), 0:(rsz - 1));
end

% inverse transformation
if opts.transinv
    dmrc.InverseTrans = inv(trf);
end

% set origin
if opts.origin
    dmrc.Origin = css;
end

% set resolution
dmrc.ResX = rsx;
dmrc.ResY = rsy;
dmrc.ResZ = rsz;

% set trans
if opts.trans
    dmrc.Trans = trf;
end
